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' In this paper we compare numerical resuhs for the ground state of the 



Hubbard model obtained by Quantum-Monte-Carlo simulations with results 
from exact and stochastic diagonalizations. We find good agreement for the 



' ground state energy and superconducting correlations for both, the repulsive 



and attractive Hubbard model. Special emphasis lies on the superconducting 
correlations in the repulsive Hubbard model, where the small magnitude of 
' the values obtained by Monte-Carlo simulations gives rise to the question, 

whether these results might be caused by fluctuations or systematic errors 
of the method. Although we notice that the Quantum-Monte-Carlo method 
has convergence problems for large interactions, coinciding with a minus sign 
problem, we confirm the results of the diagonalization techniques for small 
and moderate interaction strengths. Additionally we investigate the numer- 
ical stability and the convergence of the Quantum-Monte-Carlo method in 
the attractive case, to study the influence of the minus sign problem on con- 
vergence. Also here in the absence of a minus sign problem we encounter 
convergence problems for strong interactions. 



I. INTRODUCTION 



The Quantum-Monte-Carlo- Algorithms (QMC) are a well established tool in theoretical 
many particle physics 0. They have been successfully used to investigate the spin 
and charge degrees of freedom for the attractive and repulsive Hubbard model |0. There 
are several kinds of QMC algorithms. We focus here on the Projector Quantum-Monte- 
Carlo- Algorithm (PQMC) 0]. In the case of the repulsive Hubbard model (RHM) the sign 
problem severely restricts the parameter space one can calculate. One goal of this paper 
is to show, what parameters of the RHM can be calculated using state of the art QMC- 
methods and which of the observables still give meaningful results when a sign problem 
occurs. Concerning the minus sign free attractive Hubbard model ( AHM ) we also want to 
distinguish numerical instabilities and convergence problems from well-" behaved" runs. 

Inspired by results of de Raedt et al. for the Heisenberg model we additionally inves- 
tigate whether the choice of the elementary move in the Monte-Carlo procedure will have 
influence on the results. 

We show that exact diagonalization (ED) 0, the stochastic diagonalization(SD) [0 and 
PQMC 1^ give consistent results for a broad range of the simulated parameters. 

In our paper we investigated a modified version of the Hubbard model the t-tp-Hubbard 
model that is defined on a square lattice by the Hamiltonian: 

H=-t + S^^io) (4aCja + c]aCia) + II ^i^^il (l) 

The number operator rii^ and the fermion creation operator cj^ are defined as the ones with 
spin a at site i. a is the summation of spin up and spin down particles. The summation 
with respect to j) is over the nearest neighbor pairs, whereas the summation with respect 
to is over the next nearest neighboring pairs. Taking tp = reproduces the usual 

single band Hubbard model [Q. The hopping parameter t = 1 is taken as the energy unit 
throughout this paper. 
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II. NUMERICAL METHODS 



A. Exact Diagonalization 

A basic technique to determine the ground state properties is the exact diagonahzation. 
Here one writes the Quantum problem genproblem 

iy|^) = ^|^) (2) 

and searches for the lowest eigenvalue and eigenvector ( eigenstate ) of the Hamilton-operator 
H. We use for this paper the well known Lanczos-method 0. 

As basis states for the matrix-formalism of the Hamiltonian we use plane waves, i. e. we 
take basis-states of the form: 

•= 4it42t • • • 4„T4ii4^i • • • 4;ilo)' (3) 

where |0) is the vacuum state. To reduce the amount of the basis states we apply the 
translation symmetry, and restrict the calculation to the subspace where the total momentum 
K := Yl,{ki + k'j) = 0. As initial-vector for the Lanczos-procedure we use a random vector 

i 

with the norm 1. 

As Lanczos iteration we use 

Hqj = Pj^iqj-i + ajQj + PjQj+i with PoQo = (4) 

with Qj being the Lanczos-vector in the j-th iteration. 

The limitation of the exact diagonalization is the amount of memory that is dramatically 
increasing with increasing system-size. For a lattice with N = Ux-Uy points and n^. particles 
with spin up and rif. particles with spin down the Hilbert-space has the dimension D: 

This means, that a system with = 4 ■ 4 = 16 lattice-points and rig = 5 electrons per spin 
has the dimension D = {^^^ = 4368^ ~ 1.9- 10'' and can be solved with the Lanczos-method 
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on state of the art computers. But a system with = 6 ■ 6 = 36 and Ue = 13 electrons per 
spin direction has D = (^^fj ^ (2.3108 ■ 10^)^ ^ 5.3398 ■ 10^^ and can not be solved by exact 
diagonalization with computer systems avaible today. 

One method that circumvents some of these limitations is the stochastic diagonalization 
method. 

B. Stochastic Diagonalization 

For the stochastic diagonalization one also writes the physical quantum problem as an 
eigenproblem to extract the ground state properties. But one does not use all basis-states 
0. From the view of memory consumption it would be possible to handle all system-sizes 
with this technique. But due to a dramatic growth of the CPU-consumption with system- 
size it is so far only possible to determine the Hubbard-model for systems up to 8 x 8 sites 
for small and intermediate interactions. 

Sketching the basic ideas of the algorithm we consider n basis-states, which are or- 
thonormal. The matrix ( Hamiltonian ) for these n basis-states is then denoted as H^"'\ 
The problem is to find the eigenvalues -E^"'*. It is assumed, that they are ordered: 
^(n) < ^(n) < ... < ^(n). They can be obtained from H^"^^ with an unitary transformation 

^(n) ^ Jjin) -fjin)jj{nf ^ (^g>) 

where i?*-"-' is a diagonal matrix with the eigenvalues Ei"'\ • • • ? En{n). There are several 
standard algorithms for obtaining this unitary transformation. In the classical Jacobi algo- 
rithm 1^ for example the transformation ?7^"^ is a product of plane rotations U^"'''^\im, jm)- 

lim l[U^^'^\i^^,jJ. (7) 

m 

Each plane rotation annihilates two off-diagonal elements H'f^'"^^ and H^^'"^^ , but also zeros 
in the off-diagonal elements can become nonzero after a rotation. Only the sum over all 



4 



off-diagonal elements converges to zero. This means that, it is very important for the 
convergence, in which order the rotations are performed. Usually one always annihilates the 
largest off-diagonal element in the following step. It is impossible to perform the infinite 
number of plane rotations in equation |^ with a computer, which is necessary for the exact 
diagonalization of H^'^\ Therefore one only performs so many plane rotations, until the 
absolute value of the largest off-diagonal matrix elements is smaller than a threshold Tr. 

In the modified Jacobi scheme one only computes the smallest eigenvalue of i. e. 
one keeps im to be 1. This has the advantage, that it is not necessary to store the whole 
matrix and to make less rotations. One can prove [0 that the diagonal-element h[^i is the 
smallest eigenvalue of the matrix H^'^\ if the rotations are performed as described in 0. 

A short description of the algorithm is given by: 

1. choose one basis state as an initialization. 

2. generate a set of new trial-states. 

3. search the best trial-state for expanding the [n x 'n,)-matrix to an {{n -|- 1) x (ra -|- 1))- 
matrix. The best trial-state, is the state, for which the energy reduction Am with the 
unitary transformation U^^^^'^\l,n + 1) is the largest. 

4. if the reduction Am of the best state is bigger than a threshold Ta, than this best state 
is added to the matrix and we apply plane rotations f/('^+i'™) ( m > 1 ) to the matrix 

until the size of all off-diagonal elements f/j""^^'™') with j = 2, . . . , (n -|- 1) is 
smaller than a given threshold Tr. 

5. if the reduction A^, of the best trial-state is smaller than the threshold Ta reduce the 
threshold Ta- 

6. if the convergence criteria are satisfied stop the iterations. 

7. goto step 2. 
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One problem of the algorithm is: what are good convergence criteria? One is the conver- 
gence of the ground state energy e[^^ ( the smallest eigenvalue of the matrices ) with 
increasing number of states n. Others are the convergence of the physical properties one is 
interested in, like correlation-functions. An additional good criterion for the convergence of 
the algorithm is the norm Ni^'"^^ of the first column of the matrix H^^'^^i: 



\ 



n „ 



(8) 



C. The Projector Quantum Monte Carlo Method 

The QMC-algorithm used in this paper to investigate the ground state properties of the 
Hubbard model is the Projector-Quantum-Monte-Carlo (PQMC) algorithm. The basic idea 
of this algorithm is to use a projection operator e~^^ to extract the ground state |\t'o) of the 
Hamiltonian H from the trial wave function |\E't) H]: 



l^o) = lim c-^^I^t) (9) 
The trial wave function \^t) is chosen to be a Slater determinant. 



Sorella et al. [jTO[ have developed a stable algorithm for this purpose. We give a short 



description to introduce the notation. The Hubbard Hamiltonian (|l]) then reads as : 

H = Hkin + Hhub (10) 

Hhub = UY^ ni-^nii 

i 

Using the Trotter-Suzuki formula the exponential of the Hubbard Hamiltonian is 
rewritten as: 

m— >oo V / ^ ' 

For the interaction term Hhuh we introduce the discrete Hubbard Stratonovich transfor- 
mation |]ll| to eliminate the quartic term as: 
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(T=±l 



(12) 
(13) 
(14) 



<T=±1 



where the following definitions have been used: 

nil + ^ii ~ 1 if U < 
+ nil if [/ > 



n,; = < 



and rhi = < 



nil + nii-l 



nil - 



if < 
if [/ > 



(15) 



and 



cosh(A) = eil^l 
_ ^ 
m 



(16) 
(17) 



The elements of the discrete Hubbard Stratonovich field cr can only take the values ctj = ±1. 
The Hubbard Stratonovich field is often referred to as Ising-spins or cr-spins. 
The projection of the ground state then reads as: 



l=l,m 

cr 

The expectation value of an observable A in the ground state is now given by: 



(18) 



(19) 



Normalizing the expectation value ((T|A|cr') with {a\a') and using the following definitions 



w{a, a') = 





A 









{(T\(T 



E|(<7k')| 



(20) 
(21) 



the expectation value reads as: 



Y.w{a, (7')sign(((7|(7'))(A) 

m)-- 7 ^ ■ (22) 

E w(cr, cr')sign (cr|cr')) 

The sum over cr, cr' is now calculated by an importance sampling Monte-Carlo procedure. 
The algorithm can be sketched in the following way: 

1. create a new configuration (cr, cr')„+i from the old configuration (a, a')n by an elemen- 
tary move. 

2. calculate the transition probability 

3. calculate a random number z between and 1 

4. if 2; < P^(cr, cr')„ ^ (cr, cr%+i j accept the new configuration (cr, cr%+i else keep the 
old configuration {(T,(7')n 

5. with this procedure one generates the Monte Carlo configurations with the probability 
'u;(cr, cr'). The expectation value of an observable A is then given by: 

Esign(w(a, a')){A) 

«-^» = '-'—^) (^^) 

with the average sign (sign): 

(sign) = ^ sign(w(cr, cr')) (25) 

acr' 

In order to use the determinant {(t\(t') as a probability one has to take the absolute value. 
The sign of this determinant then becomes part of the expectation value which causes the 
so called Minus-Sign problem. 

We implemented three different kinds of elementary moves, acting on the N x m array 
of the (j-configurations (a-spins), where N is the system-size and m the number of Trotter 
slices. We are motivated by results obtained for QMC-Simulations of the Spin | Heisenberg 

8 



model where evidence was presented that the elementary move influences the convergence 
of the algorithms and can even lead to wrong results . 

1. The Cluster move flips a certain number S of randomly chosen a-spins. The criteria 
for the choice of the number S was, that the acceptance rate of these moves was at 
least 40%. 

2. The Slice move flips a certain number 5* of randomly chosen a-spins only on a single 
Trotter slice. By visiting the slices in consecutive order and storing intermediate results 
of the calculation the amount of computations can be reduced significantly. Details of 



the implementation will be given elsewhere [|T2 



3. The Single Spin elementary move only flips one randomly chosen single a-spin. This 
allows a signiflcant simplification of the algorithm [|13] . 



For the efficient implementation of the algorithm the elementary move is more than a 
minor detail and the execution time depends significantly on the choice of the elementary 
move P^ . 

QMC simulations of the Hubbard model are in general very time consuming. The calcu- 
lations were carried out on a parallel computer. Details of the parallelized algorithm have 
already been published elsewhere [jl^ . 

Finally we define two important quantities, the Monte-Carlo steps and the error bars in 
the QMC algorithm. One Monte Carlo step (MCS) is completed when the program has tried 
to flip all iV ■ m a-spins ( N is the size of the lattice and m is the number of Trotterslices ) . 

If after one run through the lattice the acceptance rate is too low, additional runs are 
performed until the acceptance rate is reached. Second the error bars are taken over bins, i.e. 
typically 50 measurements are sampled in one average A. The error bar is now calculated 
as the fluctuations of these averages A. 
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III. OBSERVABLES 



As a check whether the different programs are correct we compare different observables. 
The definitions of these observables will be presented in this section. 

The observable commonly used for the check of a correct ground state determined with 
a method is the ground state energy per lattice site Eq/N. 

A common way to study superconductivity in the Hubbard model is to examine the 
two particle correlation functions for the occurrence of Off Diagonal Long Range Order 
( ODLRO ) in the model |]r3| . In order to do this we calculate the full correlation function 
with dx2-y2 symmetry [[T^ 



that measures superconducting correlations as a function of lattice vector r. The phase 
factors gs,gs/ are 1 in x-direction and —1 in y-direction. The sum over i goes over the 
whole lattice. The sums with respect to 6 and 6' are the independent sums over the nearest 
neighbors of i. 

When studying the superconducting properties of small Hubbard-Clusters the full cor- 
relation function is not a very appropriate measure as this expectation value also contains 
contributions from the one-particle Greens functions C^'-"'^^^{r) = ^ J^iiclaCi+ra) ■ princi- 
ple these one particle contributions can be neglected as they decrease to zero with increasing 
distance \r\. But when studying small clusters their influence on the results has to be taken 
into account ||16|. Therefore we study the vertex correlation function 



C^mXr) = Cjf , (r) - 5: {gs9'sC-''%r)C-''%r + 5'- 5)) (27) 

^ a .a ^^^^ 



i,d,d' 
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IV. COMPARISON 



In this section we compare the results obtained with ED, SD and PQMC We restrict 
ourselves to the case of the 4x4 system. This is the largest square lattice that can be 
solved with the exact diagonalization. The goal is to check the correct convergence of the 
QMC-algorithms. By this we mean for the PQMC, that we check, whether the Monte Carlo 
procedure converges correctly and whether the projection parameter 6 and the Trotter de- 
composition m are chosen sufficiently large to obtain the true ground state. Special emphasis 
lies on various technical problems that come with the application of the QMC method as 
there are the minus sign problem, numerical instabilities and statistical fluctuations. In the 
case of the SD we investigate whether the number of basis states accumulated is sufficient 
to be a good approximation of the "true" ground state. 



Inspired by the results for the RHM |jT^ for larger lattice sizes we first turn our attention 
to the RHM. The most common indicator for reaching the ground state is the ground state 
energy per lattice site Eq/N. In table | and |I| we present a comparison of the ground state 
energies of the ED, SD and the three different implementations of the PQMC-algorithm 
(section II. C). The error bars on the QMC results can be estimated to be 0.25%. The 
result of ED, SD and QMC are the same within these error bounds. We find excellent 
agreement between these different methods for the ground state energy even when a minus 
sign problem is present in the calculations. This can be seen from table |Tl| where we list the 
average sign ( (sign), eq. |2^) for the simulations presented in table | and |l[ The (sign) has 
the same value in all three PQMC-algorithms we used. The SD method gives in contrast 
to the PQMC algorithm an upper bound for the exact ground state energy Eq. Therefore 
comparing two energies obtained with SD it is clear, that the lower energy is closer to the 
exact energy. For the PQMC the error bar gives a bound for the ground state energy, but 
it is not possible to deduce from the PQMC value whether the energy is above or below the 
exact energy. 

It has turned out that comparing ground state energies is not a very sensitive indicator 
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whether the PQMC has reached the true ground state |jT8[, |[T9[ . 



Our interest hes on the expectation values of the superconducting correlation function 
for the RHM in the (ia,2_y2-channel. 

In figure 01 we present the comparison of Cj"^' (r) (eq. pUf) for five different methods 
namely the ED, SD, single-PQMC, slice-PQMC and cluster-PQMC. For a 4 x 4 Hubbard 
system with tp = —0.22, A^e = 5 + 5 and [/ = 2, 4, 5 and 6 respectively we observe perfect 
agreement of all five programs. The points for all five methods fall on identical positions 
and cannot be distinguished. 

As we outlined in the previous section the more appropriate indicator for the existence 
of ODLRO in small clusters is the vertex correlation function Cj'^^^'^^{r) (eq. [27| ). The 
absolute values of this vertex correlation function compared to the full correlation function 
turned out to be very small for the RHM giving rise to the question whether these 
small results might be artifacts of the various problems that come with the application of 
the QMC-method. 

In figure ^ a-d we present a systematic investigation of the 4 x 4-system with 5 + 5 
electrons for various values of U. In contrast to the case of the full correlation function 
Cj"^' (r) we observe that in the case of the vertex correlation function CY'^J'^'^^ (eq. PTI) only 
simulations with U = 2 and U = 4 show again a perfect agreement. For larger interactions 
U = 5 and f/ = 6 we find strong fluctuations of the different PQMC methods coinciding 
with the occurrence of a minus sign problem in the simulations (table |ITTD . Thus we have 
shown that whereas the ground state energy and the full correlation function seem perfectly 
converged this is no proof that the vertex correlation function is converged in these runs. 
The fluctuations in the correlation functions have increased dramatically with interaction U 
despite the significant higher number of MCS, as shown in table |V[ 

We now turn our attention to the AHM. Here we consider both the cluster and the single 
spin dynamics in the PQMC algorithm. The AHM shows superconductivity in the on-site 
s-wave channel. In this case there is no summation over the nearest neighbors, i.e. no sum 



with respect to 5 and S' in the eq. |2g and Of course the phasefactors gs = gs' are equal 
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1. 

First we compare the ground state energy for PQMC with ED and SD (table |V]). Note 
that for the SD deviations from the ED values develop with increasing interaction strength 
\U\. For interactions U = —6 and U = —8 the SD algorithm has convergence problems. For 
these high values of \U\ the k-space representation of the basis states is no longer appropriate. 
This can be seen from figure |^. The SD algorithm is able to find almost every important 
state (i.e. all states with a large weight) for the interaction U = —2. In the interval 0.0001 
to 0.00032 there are 6276 states in the ground state of the exact diagonalization. 6012 states 
are in the same interval in the basis of important states of the SD. For U = —2 the SD has 
collected almost every state with a weight larger than 10~^. 

The number of states with a relatively large weight increases for stronger interactions 
\U\. In the interval 0.001 to 0.0032 the number of states in the ground state of the exact 
diagonalization goes from 1 082 {U = —2) to 24 003 {U = —8) (fig. Therefore more states 
are necessary as a good approximation of the ground state in the SD. These convergence 
problems can be demonstrated even clearer in table |V1[ This table shows additionally to 
the number of states Nsd in the basis of the SD (plotted in fig. |^) the norm 
the first column (eq. ^) and the threshold for the acceptance of new states. Because 
the CPU-time is limited [l^, [l^ in practice the number of states Nsd can not exceed to 



some hundred thousand. From table |VI| one can see, that Nsd and A?^}"'"*^ as well as are 
increasing with \U\. For practical purposes the SD is not able to find a good approximation 
of the exact ground state beyond a certain interaction. We were not able to resolve these 



problems by the use of a real space representation of the basis states pO |. 

As a second step we compare the superconducting correlation function with on-site s- 
wave symmetry for the AHM. Figure ^ shows C^"" and |^ shows C^f^^'^^{r) for U = —2, —4, 
—6 and —8 for a 4 x 4 system with tp = 0. The comparison of the vertex correlation function 
with the Lanczos method clearly indicates that for U = —6 and especially for U = —8 the 
cluster program shows deviations from the ED values. These deviations are comparable to 
the case of the RHM, where we also observe convergence problems for increasing interaction 
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U. We tried to find out where these convergence problems might come from. On the one 
hand we investigated the dependence on the projection paramter 6. In figure |^ and | we 
show the dependence of the cumulated plateau of the vertex correlation function C^^^'*'^^{r) 
and of the ground state energy per site Eq/N of the projection parameter O. It clearly 
indicates that a projection parameter of = 8 as it was used throughout the rest of the 
paper is sufficient to reach the ground state. Note that especially the ground state energy 
does not depend very strongly on the projection parameter. 

Furthermore we investigated whether numerical instabilities might cause the convergence 



problems. In table |VII| we compare the ground state energy for different numbers of stabi- 
lizations nstab. To reduce numerical instabilities the product of m imaginary time slices is 
being stabilized by a modified Gram-Schmidt orthogonalization . 



A stabilization is performed every nstab slices of the Trotter decomposition. The results 
agree perfectly so we conclude that the simulations were numerically stable. 

Finally we checked whether the MC process depends on the starting seed of the random 
number generator. Figure ^ shows the ground state energy for different seeds. Here we also 
cannot discover any unusual behaviour. 

Thus even in the absence of a minus sign problem for the AHM the PQMC algorithms 



have convergence problems. An even more striking result is seen comparing table |3 and 



VIII| . The number of MCS needed for convergence both in the RHM as well as in the AHM 
increase with interaction strength \U\ for both the single and the cluster algorithm in the 
same manner. As a comparison we show in fig. ^ the results for unconverged runs. The 
corresponding number of MCS (PQMC) and number of basis state (SD) are listed in table 

E 

Given only limited CPU power, we come to the result, that even when there is no minus 
sign problem as for the example in the case of the PQMC algorithm for the attractive 
Hubbard model or of the SD the presented numerical algorithms cannot handle values of U 
beyond a critical interaction, because of convergence problems. 



14 



V. CONCLUSION 



We presented a detailed comparison of the results of exact diagonalization, stochastic 
diagonalization and the PQMC. For the PQMC we implemented three different elementary 
moves, namely the single spin flip move, the slice move and the cluster move. We showed 
that there is perfect agreement between all five methods for the ground state energies and 
the full correlation function for the RHM and the AHM. Also for the vertex correlation 
functions, that are very small in magnitude in the repulsive case, the exact diagonalization 
results can be confirmed by PQMC and stochastic diagonahzation. But for large interactions 
it is necessary to perform dramatically more MCS in the PQMC respectively to collect much 
more states in the SD. The observation that with an increase of |[/| the amount of CPU-time 
is enhanced dramatically leads to the problem that these calculations cannot be performed 
with state of the art computers for large interaction strengths \U\. This increase in CPU-time 
was found both, in the RHM where it was expected due to the occurrence of the well-known 
minus sign problem and in the AHM that is minus sign free. 
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TABLES 



u 


Eed/N 


Esd/N 


E single/ N 


E slice /N 


Eduster 


2.0 


-1.336059 


-1.3360 


-1.337 ±0.001 


-1.335 ± 0.001 


-1.332 ± 0.003 


4.0 


-1.223808 


-1.2237 


-1.224 ±0.002 


-1.224 ± 0.002 


-1.22 ±0.01 


5.0 


-1.182000 


-1.1819 


-1.180 ±0.003 


-1.184 ± 0.003 


-1.19 ±0.02 



TABLE I. Comparison of the ground state energies per lattice site for a Hubbard model of 
system size A = 4 x 4, number of electrons A^e = 5 ± 5, diagonal hopping-matrix element tp = 0.0 
for various interactions U. In column 1 we show the interaction U, in column 2 we present the 
results for the exact diagonalization (Eed/N), in column 3 we show the stochastic diagonalization 
results (Esd/N), in column 4-6 we present the results obtained with the PQMC- Algorithms for 
different elementary moves. Details of these moves are given in the text. 



u 


Eed/N 


Esd/N 


Esingle/N 


E slice /N 


Eduster /N 


2.0 


-1.230034 


-1.2300 


-1.2301 ± 0.0003 


-1.2306 ± 0.0008 


-1.229 ±0.001 


4.0 


-1.126160 


-1.1261 


-1.125 ±0.003 


-1.125 ±0.002 


-1.13 ±0.02 


5.0 


-1.088907 


-1.0887 


-1.089 ±0.004 


-1.087 ±0.004 


-1.10 ±0.02 


6.0 


-1.058717 


-1.0581 


-1.061 ±0.005 


-1.058 ±0.007 


-1.06 ±0.02 



TABLE II. Comparison of the ground state energies per lattice site for a Hubbard model of 
system size A = 4x4, number of electrons Ag = 5±5, diagonal hopping-matrix element tp = —0.22 
for various interaction U. In column 1 we show the interactions U, in column 2 we present the 
results for the exact diagonalization (Eed/N), in column 3 we show the stochastic diagonalization 
result (Esd/N), in column 4-6 we present the results obtained with the PQMC- Algorithms for 
different elementary moves. Details of these moves are given in the text. 
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U = 2 


U = A 


U = 5 


U = 6 


tp = -0.22 


1.0 


0.92 


0.69 


0.42 


tp = 0.0 


1.0 


0.99 


0.95 


0.85 



TABLE III. Average sign (sign) for 4 x 4 Hubbard cluster with ATg = 5 + 5, (sign) is in the 
given precision the same in all three implementations of the PQMC algorithm. 



U cluster single slice 

2.0 790 000 49 000 23 400 

4.0 890 000 79 000 50 000 

5.0 4490 000 199 000 65 200 

6.0 5 390 000 599 000 72 900 



TABLE IV. Monte Carlo Sweeps MCS for cluster and single spin flip algorithm for the a 4 x 4 
Hubbard model with Ng = 5 + 5, tp = —.22 and interaction U of the runs in figure 1, 2, 4 and 5 



u 


Eed/N 


Esd/N 


E single/ 


E duster /-^ 


-2.0 


-1.731689 


-1.7316 


-1.731 ± 0.003 


-1.732 ± 0.006 


-4.0 


-2.045849 


-2.0453 


-2.045 ± 0.002 


-2.048 ± 0.01 


-6.0 


-2.458782 


-2.4568 


-2.460 ± 0.004 


-2.47 ± 0.01 


-8.0 


-2.956890 


-2.9545 


-2.952 ±0.01 


-2.969 ± 0.05 



TABLE V. Comparison of the ground state energies per lattice site for an attractive Hubbard 
model of system size AT = 4 x 4, number of electrons ATg = 5 -|- 5, diagonal hopping matrix element 
tp = Q for various interactions U. In column 1 we show the interaction U, in column 2 we present the 
results for the exact diagonalization {Eed/N), in column 3 we show the stochastic diagonalization 
result {Esd/N), in column 4-5 we present the results obtained with the PQMC- Algorithms for 
different elementary moves. Details of these moves are given in the text. 
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u 




NsD 




Ta 


6.0 


-.22 


307 324 


0.0491 


0.381 • 10-"^ 


5.0 


-.22 


285 872 


0.0426 


0.381 • 10-^ 


4.0 


-.22 


228 172 


0.0334 


0.381 • 10-^ 


2.0 


-.22 


141 742 


0.0062 


0.954 • 10"^ 


1.0 


-.22 


60 638 


0.0020 


0.238 • 10"^ 


-1.0 


0.0 


20 974 


0.0040 


0.954 • 10-^ 


-2.0 


0.0 


28 034 


0.0094 


0.298 • 10"^ 


-2.0 


0.0 


115 304 


0.0140 


0.191 • 10"^ 


-4.0 


0.0 


179 755 


0.1117 


0.153 • 10-6 


-6.0 


0.0 


253 275 


0.1297 


0.153 • 10-6 


-8.0 


0.0 


263 862 


0.1514 


0.153 • 10-6 



TABLE VI. Number of states NgD used for the stochastic diagonahzation method in a 4 x 4 
system with iVg = 5 + 5 electrons and = for different interaction U . N^'^^ is the norm of the 
first column in the transformed matrix H^"'\ and Ta is the threshold for the acceptance of a new 
basis states to the basis of important states. 



nstab MCS Eduster 

1 1745000 -2.45282 ± 0.00607 

2 1745000 -2.45282 ± 0.00607 
4 1745000 -2.45282 ± 0.00607 
8 1745000 -2.45282 ± 0.00607 

16 1745000 -2.45282 ± 0.00607 

TABLE VII. Groundstate energy for 4 x 4 Hubbard cluster with = 5 + 5, tp = 0.0 and 
U = —6 for different numbers of stabilizations nstab. 
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u 



cluster 



single 



-2.0 39 000 4000 

-4.0 189 000 9 000 

-6.0 2 909 000 109 000 

-8.0 7 750 000 149 000 

TABLE VIII. Monte Carlo Sweeps MCS for the cluster and the single spin flip algorithm for 

a 4 X 4 Hubbard model with Ne = 5 + 5, tp = 0. The interactions U are the same of the runs in 
figure 4 and 5. 



U Eed/N Esd/N Esingle/N Eauster/N 

-2.0 -1.731689 -1.7314 -1.731 ±0.003 -1.732 ±0.006 

-4.0 -2.045849 -2.0419 -2.035 ±0.003 -2.078 ± 0. 

-6.0 -2.458782 -2.4451 -2.457 ± 0.002 -2.439 ± 0. 

-8.0 -2.956890 -2.9432 -2.952 ± 0.001 -2.977 ± 0. 

TABLE IX. Comparison of the ground state energies per lattice site for an attractive Hubbard 
model of system size AT = 4 x 4, number of electrons Ag = 5 + 5, diagonal hopping matrix element 
tp = for various interactions U. Here we compare unconvcrgcd runs. In column 1 we show 
the interaction U, in column 2 we present the results for the exact diagonalization {Eed/N), in 
column 3 we show the stochastic diagonalization result (Esd/N), in column 4-5 we present the 
results obtained with the PQMC- Algorithms for different elementary moves. Details of these moves 
are given in the text. 
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U state in SD cluster single 

-2.0 12122 39 000 4 000 

-4.0 23 510 50 000 5 000 

-6.0 32 587 1000 000 49 000 

-8.0 46 583 1500 000 109 000 



TABLE X. Number of basis states in the SD and Monte Carlo Sweeps MCS for the cluster 
and the single spin flip algorithm for a 4 x 4 Hubbard model with Ng = 5 + 5, tp = and attractive 
interaction U of the unconverged runs of figure 6. 
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FIGURES 

FIG. 1. Comparison of the full correlation function with d^2_y2 symmetry Cj"'' (r) for a 4 x 4 
Hubbard cluster with iV^ = 5 + 5, ip = -0.22 and [/ = 2, 4, 5 and 6. ( 6 = 8, r = 0.125 ) 

FIG. 2. Comparison of the vertex correlation function with d^2_y2 symmetry (jj<^^t<^^(^f'^ for a 
Ax 4 Hubbard cluster with Ng = 5 + 5, tp = -0.22 and [/ = 2, 4, 5 and 6. ( 6 = 8, r = 0.125 ) 

FIG. 3. Comparison of the weights of the basis states in the exact diagonalization (ED) and 
the stochastic diagonalization (SD) in a 4 x 4 system with Ne = 5 + 5 electrons for the interaction 
?7 = -2 and J7 = -8 and tp = 0.0. 

FIG. 4. Comparison of the full correlation function with on-site s-wave symmetry Cg^^\r) for 
a 4 X 4 Hubbard cluster with A^e = 5 + 5, = and [/ = -2,-4, -6 and -8. ( = 8, r = 0.125 ) 

FIG. 5. Comparison of the vertex correlation function with on-site s-wave symmetry C^^''*^^(r) 
for a 4x4 Hubbard cluster with = 5+5, tp = and i7 = -2, -4, -6 and -8. ( G = 8, r = 0.125 ) 

FIG. 6. Comparison of the vertex correlation function with on-site s-wave symmetry C^f^^^^{r) 
for a 4 X 4 Hubbard cluster with Ne = 5 + 5, tp = and U = -2, -4, -6 and -8. ( 6 = 8, 
r = 0.125), The runs are unconverged meaning that for \U\ > —2 the number of Monte Carlo 
Steps was chosen significantly lower than in figure 5. Details are given in table X 

FIG. 7. Investigation of the behaviour of the vertex correlation function with on-site s-wave 
symmetry CY^J'^^^{r) for a 4 x 4 Hubbard cluster with A^e = 5 -|- 5, tp = and —6 for 9 = 1,2,4,8 
and 16. (r = 0.125) Dashed line is the value of the exact diagonalization. 

FIG. 8. Investigation of the behaviour of the ground state energy for a 4 x 4 Hubbard cluster 
with iVe = 5 -I- 5, tp = and -6 for = 1,2,4,8 and 16. (r = 0.125) Dashed line is the value of the 
exact diagonalization. 



21 



FIG. 9. Investigation of the behaviour of the ground state energy for a 4 x 4 Hubbard cluster 
with ATg = 5 + 5, = and U = —6 for different seeds of the random number generator. (0 = 8, 
r = 0.125) Dashed line is the value of the exact diagonalization and the full line is to guide the 
eyes. 
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